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Semidiscrete Galerkin Modelling of Compressible Viscous Flow 
Past a Circular Cone at Incidence 

By 


Andrew J. Meade Jr. 


Mechanical Engineering 

CtWL ^/lT 

Chairman, Thesis Committee 


Abstract 


A numerical study of the laminar and compressible boundary layer, about a circular 
cone in a supersonic free stream, is presented. It is thought that if accurate and efficient 
numerical schemes can be produced to solve the boundary layer equations, they can be 
joined to numerical codes that solve the inviscid outer flow. The combination of these 
numerical codes is competitive with the accurate, but computationally expensive, Navier- 
Stokes schemes. 

The primary goal of this study is to develop a finite element method for the calcula- 
tion of three dimensional compressible laminar boundary layer about a yawed cone. The 
proposed method can, in principle, be extended to apply to the three dimensional boun- 
dary layer of pointed bodies of arbitrary cross section. 

In the present thesis the three dimensional boundary layer equations governing super- 
sonic free stream flow about a cone are examined. The three dimensional partial differen- 
tial equations are reduced to two dimensional integral equations by applying the Howarth, 
Mangier, Crocco transformations, a linear relation between viscosity, and a Blasius-type of 
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similarity variable. This is equivalent to a Dorodnitsyn-type formulation. The reduced 
equations are independent of density and curvature effects, and resemble the weak form of 
the two dimensional incompressible boundary layer equations in Cartesian coordinates. In 
addition the coordinate normal to the wall has been stretched, which reduces the gradients 
across the layer and provides high resolution near the surface. 

Utilizing the parabolic nature of the boundary layer equations, a finite element 
method is applied to the Dorodnitsyn formulation. The formulation is presented in a 
Petrov-Galerkin finite element form and discretized across the layer using linear interpola- 
tion functions. Linear functions were chosen for ease of programming while also providing 
adequate accuracy. The finite element discretization yields a system of ordinary differential 
equations in the circumferential direction. The circumferential derivatives are solved by an 
implicit and noniterative finite difference marching scheme. 

Solutions are presented for a 15° half angle cone at angles of attack of 5° and 10°. 
The numerical solutions assume a laminar boundary layer with free stream Mach number 
of 7. Results include circumferential distribution of skin friction and surface heat transfer, 
and cross flow velocity distributions across the layer. 
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1. INTRODUCTION 


The supersonic flow over sharp conical bodies at angles of attack has been a topic of 
considerable interest, both theoretically and experimentally, for the past three decades. This 
interest stems from the following properties. 

(i) The geometry of conical bodies possess a high lift to drag ratio at hypersonic speeds 
and excellent trajectory accuracy. This makes conical bodies prime re-entry vehicles. 

(ii) There is a requirement in future aircraft and missiles, which possess conical forebo- 
dies, for greater control at high angle of attack flight. The characteristics of high angle 
of attack flight are strong viscous-inviscid interaction and three-dimensional flow 

separation. 

(iii) The theoretical assumption of supersonic conical flow allows for simplification in 
analysis while it also reveals detailed flow behavior that is typical of more complex 

geometries. 

The theoretical investigations of the flowfield can be divided into two main classes, 
corresponding to the inviscid and viscous models. The viscous models require the solution 
of an approximation of the Navier-Stokes equations for the simulation of the flowfield. Hel- 
li well and Lubard [1] solved the approximation of the full equations (ignoring streamwise 
viscous diffusion) for laminar flow about a cone, at small to moderate angles of attack. 
McRae and Hussaini [2] also used the same Navier-Stokes approximations to solve for both 
laminar and turbulent flow, at moderate incidence. Degani [3] has solved a similar set of 
equations for the cone at high angles of attack for turbulent flows. All of the solutions men- 
tioned model most of the physical mechanisms and should give a more accurate prediction 
of the flow field than the inviscid models. However, the greatest drawback these viscous 
models possess, are their requirements for large amounts of computer time and storage. It 
is found that the inviscid models, because of their relative speed and simplicity, are still in 


wide use in theoretical investigations. 
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For inviscid analysis the flow field is considered conical. While the flow is in a real 
ihrf' ('-dimensional, all flow variables depend on only two coordinates 6 and (Fig. 
la and b). Since the partial differential equations are nonlinear, with unknown boundary 
conditions at the outer shock, a numerical method must be used in order to obtain solutions. 
Historically the first attempt to obtain a numerical solution for the flow was made by Stone 
[4J.|5], Many of the inviscid models used for conical bodies make use of vortex filaments 
[6], point vortices [7], or the popular Euler equations. As a less expensive alternative to the 
Navier-Stokes equations, the Euler equations offer many attractive features. In contrast to 
potential methods, the Euler equations provide the correct Rankine-Hugomot shock jump 
conditions, and correctly describe the transport of vorticity. The solution of the Euler equa- 
tions about the cone have been made by a variety of methods. Reference [8] contains a 
comprehensive bibliography of the early solution schemes. The more recent methods 
include the popular shock -capturing [9] and finite volume schemes [10]. 

AH of the inviscid models considered assume that viscous effects can be approximated 
by experiments or theoretical means. The yawed cone in supersonic flow has been studied 
experimentally by Holt and Blackie [11], Tracy [12], Rainbird [13], and Yahalom [14], 
among others. Yahalom contains an excellent bibliography of experimental studies prior to 
1971. Some of the later studies include references [15],[16], and [17], Unfortunately, 
experimental data are usually expensive and difficult to obtain for the conditions of interest. 
Moderate to high angles of attack result in flow separation from the cone, and require an 
interaction process between the inviscid and the significant viscous flow regions. Experi- 
mental data are not sufficient. The most straightforward and accurate means of accomplish- 
ing the interaction, requires an iterative procedure between the external flow and a theoreti- 
cal viscous model. 

The most familiar theoretical means of accounting for the viscous effects, is the solu- 
tion of the classical boundary layer equations. The purpose of this text is to develop an 
accurate, computationally efficient, and general numerical scheme for the solution of the 
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boundary layer equations on a cone. Ultimately, it is hoped that by linking this numerical 
method to an efficient Euler equation solver, the resulting inviscid model can be made com- 
petitive with the existing viscous models. 

An early major paper on conical boundary layers, that contains many of the ideas used 
in this text, is that by Moore [18]. Moore established the equations of compressible, three- 
dimensional boundary layer flow over a yawed circular cone. Solutions by later investiga- 
tors [19]-[23], utilizing the equations of Moore, mostly used finite difference schemes, with 
the exception of Fletcher and Holt [24], The latter applied an adaption of the method of 
integral relations (specifically the spectral method) to the governing partial differential equa- 
tions, and reduced them to ordinary differential form. The method of integral relations 
(MIR), of reference [24], is equivalent to the method of weighted residuals (MWR) 
[25], [26] found in finite element literature [27], [28]. MIR permits considerably more accu- 
rate solutions to be obtained than were previously possible. However, an increase in the 
number of coefficients did not necessarily mean an increase in accuracy. Also, their method 
gives solutions up to, but not beyond, separation. 

The numerical scheme presented in this text uses a combination of a Dorodmtsyn for- 
mulation of the boundary layer equations, and a finite difference/finite element procedure 
(semidiscrete Galerkin method) in solving the boundary layer on a yawed cone in super- 
sonic flow. 

The application of MWR, and the various transformations presented in this text, 
results in a Dorodnitsyn formulation of the boundary layer equations for the cone [29], [30]. 
This formulation is known to obtain relatively accurate solutions with only a few 
coefficients defining the dependent variables across the boundary layer. The high accuracy 
in this formulation comes, in part, from the use of the streamwise velocity component u. 
This transformation, in effect, acts as a stretched surface normal coordinate. The new 
independent variable, u, minimizes the normal gradients and provides high resolution near 
the wall. In addition the boundary layer equations are reduced to a set of ordinary 
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differential equations. 

The finite element method uses piecewise continous polynomials to model dependent 
variables and averages out the approximation errors across the elements. This averaging 
allows considerable accuracy for a small number of elements. Greater accuracy can be 
obtained by simply increasing the number of elements. Since the boundary layer equations 
for the cone are parabolic, the circumferential dependence of the variables can be modelled 
by an implicit finite difference marching scheme. The surface normal dependence of the 
variables is modelled by one-dimensional elements. This arrangement is much more compu- 
tationally efficient than that using two dimensional elements, and easier to program. In the 
past, two-dimensional boundary layer flows have been computed effectively with finite ele- 
ment formulations [31], [32], However, the previous formulations have never been applied to 
the cone problem, nor taken full advantage of the various transformations to be found in 

this text. 

The semidiscrete Galerkin method can be applied to cones with more general cross 
sections without modification. The influence of the cross section appear only through the 
external flow parameters that are found in equation (3.13). The method also possesses the 
potential capability of calculating solutions beyond flow separation. 

The present studies are restricted to attached, inviscid supersonic external flow at 
moderate angles of attack. The boundary layer is assumed laminar and the fluid is assumed 
ideal with a specific heat ratio of 1.4 (air). The only geometry considered is that for a circu- 
lar cone. 

In chapter 2 the semidiscrete Galerkin method is explained and demonstrated, using a 
one-dimensional parabolic equation. Chapter 3 describes the equations of motion used in 
the inviscid and viscous regions, and boundary and initial conditions used for the cone. In 
chapter 4 the method of weighted residuals, and the semidiscrete Galerkin method in partic- 
ular. is applied to the boundary layer equations. Chapter 5 presents the details of the 
numerical methods used to solve the equations given in chapter 4. Chapter 6 gives the 
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numerical results and includes the convergence curves, circumferential distribution of skin 
friction and surface heat transfer, and cross flow velocity distributions. Chapter 7 presents 
the summary of the text and the work in progress. 
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2. THE GALERKIN FINITE ELEMENT METHOD 


The finite element method is an approximate method of solving a wide range of boun- 
dary and initial value problems. Since there are, of course, various comprehensive refer- 
ences on the subject of finite elements [33]-[41] this chapter will only supply sufficient 
background and touch on the points needed to illustrate the method, as it will be applied to 

the cone problem. 

2.1 The Method of Weighted Residuals 

Among the more popular methods in the numerical solution of fluid dynamics prob- 
lems are finite difference, variational methods and the methods of weighted residuals [26], 
This large class of methods can be described in the following manner. Using linear opera- 
tors, it is assumed that an approximate solution to a differential equation 

L(u) = 0 (2J) 

is to be found subject to the boundary conditions 

B(u) = 0 (2 ’ 2) 

and the initial conditions 


l(u) = 0 

An approximate solution to u is introduced. 


(2.3) 


u a (x,t) = u°(x,t) + £Nj(x)aj(t) , (2 - 4) 

j=i 

where a.(t) are unknowns, N/x) are known analytic functions (trial functions), and u°(x,t) is 
chosen to satisfy all global boundary conditions (2.2) and initial conditions (2.3). Substitu- 
tion of the approximate solution into (2.1) gives an error (residual), R. To determine the 
values of a r an inner product of the residual R and a set of linearly independent weighting 
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functions w k (k = 1 N) is constructed and set equal to zero. 

Jw k Rdx = 0 (2,5) 

D 

where D is the domain of interest. 

This resembles the weak form of (2.1). As equation (2.5) is written equation (2.1) 
becomes a system of N ordinary differential equations in t. If the residual R is considered 
continuous, then as N— »°°, R should converge to zero in the mean. This implies that if the 
approximate solution u a satisfies the boundary conditions exactly then u a should converge to 
the exact solution u in the mean. In short, by setting the inner product to zero, the errors 
introduced by u a are averaged out over the domain of interest. 

The choice of the weighting function determines the various subclasses to be 
employed. The more popular subclasses are the least squares method, the method of 
moments, the collocation method, the subdomain method, and the Galerkin method. 

2.2 The Petrov-Galerkin Method 

Basically for the subclass of Petrov-Galerkin methods [42] the weighting function is 
written as 


where P k (x) is similar to the trial function Nj(x) used in u a , but with some modifications 
used to satisfy requirements on the solution. The classical Galerkin method [43] uses only 
the trial functions for weighting. Classical Galerkin methods can be thought of as a 
simplified form of the Petrov-Galerkin method. To simplify matters the Petrov-Galerkin 
method will be referred to as the Galerkin method in this text. 

To ensure that linearly independent equations exist for the solution of the values aj, 
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the trial functions N/x) and thus the weight functions must be linearly independent. 

The choice of the trial functions to be used forms the two main branches known as the 
spectral method and the Galerkin finite element method. 


2.3 The Semidiscrete Galerkin Method 

The most popular branch of the Galerkin method is the Galerkin finite element 
method. Its popularity comes from the ease in which it can be applied to nonstructural prob- 
lems. Basically the method uses low order polynomials as interpolation and weighting func- 
tions and then applies the Galerkin method within the subdomains (finite elements) formed 
from the domain of interest. 

The semidiscrete Galerkin method [28] uses finite element representation only for spa- 
tial variation. The time derivative is replaced by some finite difference operator. 

To illustrate the method and its properties, as it was applied to the cone problem, we 
will consider the one-dimensional unsteady heat conduction problem in nondimensional 
form. 


du_3^u =0 (2.7) 

at 3x 2 
for 

0 < x < 1 

where u(x,t) represents the nondimensional temperature. The initial and boundary condi- 

u(x,0) = u°(x) 
u(0,t) = 0 
u(l,t)= 1 


tions are 
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The domain of interest (0<x< 1) is divided into subdomains known as finite elements. 
As seen in figure 2 the nodal points are the element boundaries, using the variable Xj to 
mark the position of the nodes. The approximate solution for u is 

».<*.■>- <2 - 8a) 

j 

and for u° 


u a (x) = X Njiij 0 , 


(2.8b) 


where Nj is the interpolation function, and Uj and uj* are the values of u and u° respectively 
at the nodes. Using the interpolation function as the weighting function , N k , and taking 
the weak form of (2.7) we have 


' 8u a l 3 2 u a 

fN k — dx-]N k —^-dx = 0 
l di i 3x‘ 


(2.9) 


d”u 

As it stands the term — ^ cannot be represented properly. The finite element method 
9x" 

requires that the interpolation used should provide interelement continuity of derivatives of 
degree one less than the maximum that appears in the weak form of the equation of interest. 
Piecewise continuous polynomials provide continuity for only the value u a across the ele- 
ments. Also, if at all possible, we would like to employ linear interpolation functions. If 
these simple functions can be shown to be completely satisfactory, there are no compelling 
reasons to use complicated, and numerically costly, higher order functions. Linear interpola- 
tion provides continuity within the element for only a first order derivative. 

The order of the derivative can be reduced by applying the Green-Gauss theorem to 
the second integral. In one dimension the Green-Gauss theorem manifests itself as integra- 


tion by parts. Equation (2.9) becomes 
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jN 


k 


0 


du a \ dN k du a 

dx = -j - 7 — dx + N k 


at 


dx ax 


du a 

ax 


( 2 . 10 ) 


The resulting boundary term is known as the natural boundary condition. If the exam- 
ple were of Neumann type the boundary term would be replaced by the given boundary 
condition. Since the example is of Dirichlet type the natural boundary condition has been 
incorporated into the governing equation. Using the full representation of u a (2.8a) in equa- 
tion (2.10) results in 


N 

I C1 M 


J 


dUj 

dT 


N N dN; 

= -IC2 kj u j+ lN k — Uj 

J J 


( 2 . 11 ) 


where 


Cl kj =jN k N J dx 


L dN k dN; 


The semidiscrete Galerkin method has reduced the partial differential equation (2.7), in 
x and t, to a system of first order ordinary differential equtions in t. Using piecewise con- 
tinuous polynomial interpolation functions, the resulting equations will be linearly indepen- 
dent for all N. In addition, if low order polynomials are used the integrands of the equations 
will also be low order polynomials, which can be solved exactly and efficiently by low 
order Gaussian quadrature [44], Notice that the integrals are independent of t. If a uniform 
grid is used for x the values of Cl kj and C2 k j are constant. 


For linear interpolation functions at x } (Fig. 3) 


X - Xj_, 

in element A, N.= 

> x ; x j — ] 
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in element B, 


N,= 


*j+l 

Xj+l-Xj 


forx<Xj_, Nj = 0 

forx>x J+) Nj = 0 

At a specific node x k , the only nonzero contributions in the integrals come from the 
neighboring nodes, x k _] and x k+ ], and the node itself x k . Therefore Cl k j and C2 kj form tri- 
diagonal matrices. 

The boundary term becomes 


u, = 0 (2.12) 

u N = 1 

Equation (2.11) becomes 

l-l N (2.13) 

i J 

where C2 kj has been modified to absorb the natural boundary terms and renamed SC2 kj . 

The essential boundary conditions of equation (2.12) reduce the NxN arrays, Cl kj and 
SC2 kj . into (N-2)x(N-2) matrices. 

The ordinary differential equation is solved by introducing a discretized time domain 
consisting of n number of time increments, At. The time derivative is replaced by the finite 
difference form 


u J n+l -u J n _ Auj n+I 
At ~ At 


(2.14) 


where the superscript indicates the time level. The variable Auj 1 is solved using the theta 
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method [45] for integration. 


fci kl A Uj n+1 = -At(0fsC2 kj uj’ + ' + (1 - e)fsC2 kj up (2.15) 

i j J 

where 6 controls the degree of implicitness. The value of 6 =0 gives the explicit Euler for- 
ward difference scheme, while 6 = 1 gives the fully implicit backward Euler difference 
scheme. By setting 6 =0.5 the second order accurate Crank-Nicholson [45] scheme can be 
used. 

Since the governing equation is linear, it can be simply rewritten as 

£(Cl kj + 6 AtSC2 kj )A Uj n+l = -AtXSC2 kjUj n (2-16) 

j j 

For a chosen At the right hand side term is known and forms a vector of length N-2. 
The array Cl kj + 0AlSC2 kj is a (N-2)x(N-2) tridiagonal matrix that can be solved quite 
efficiently by the Thomas algorithm [46], The Thomas algorithm factorizes the tridiagonal 
matrix in O(N) operations as opposed to the 0(N 2 ) operations required by Gaussian elimi- 
nation [47]. 

To inlegrate in t, Auj^ 1 is solved for and added to the known value uj* to give uj 1 * 1 . 

The Uj n+1 is then used in the right side to solve for the next time level. Notice that by vary- 
ing the time step for required accuracy, iteration is not needed. This results in a computa- 
tionally efficient integration algorithm. 

To summarize, the application of the semidiscrete Galerkin method with the theta 
method for temporal integration, and the Thomas algorithm for matrix solution, makes for a 
very efficient computational solution of a parabolic problem. 



13 


3. GOVERNING EQUATIONS 

From a computational point of view it is more practical and convenient to classify the 
flow about the cone into two regions, the inviscid and viscous. One may then exploit the 
character of the respective regions and use efficient numerical methods in their solution. 

3.1. Inviscid Flow 

3.1.1. Equations of Motion 


In this present study the flow outside the boundary layer is assumed to be steady, 
supersonic and attached. The resulting main shock is assumed attached to the apex of the 
cone. The fluid is assumed to be inviscid and non-conducting. A spherical coordinate sys- 
tem is used with r measured from the cone apex, 6 measured from the cone axis, and <p 
measured from the windward line of symmetry (Fig. la, b). The dependent variables are 
nondimensionalized in the following manner, 



where a is the local speed of sound and a* is the critical speed of sound. The nondimen- 
sional form of the equations of motion for this flow are as follows. 


Conservation of mass: 
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JtL + L_ + — ^ — — -^2- + 2u + vcot0 = 0 (3.2a) 

30 sin0 3 p p 30 psin0 dp 


Conservation of r momentum: 


v 9u + _^L*L_ (v ! + w ! ) = 0 c- 2t » 

dO sin# dp 


Conservation of 0 momentum: 


v l^ + _^_^ + l|E+uv-w 2 cot0 =0 
dO sin# dp p dO 


(3.2c) 


Conservation of 0 momentum: 


v + + — - — & + uw + vwcot# = 0 (3. 2d) 

5# sin0 dp p sine d<p 


Conservation of energy: 


dv w dp *> dp , w dp 

v f + * — a~ v — + — — -r— 

d# sine dtp dO sine dip 


= 0 


(3.2e) 


Note that since the main shock is attached (Fig. la.b), the dependent variables are indepen- 


dent of the 7 coordinate. 


3.1.2. Boundary Conditions 

The external flow field is bounded by the body, the main shock, and the windward and 
leeward planes of symmetry ( p = 0° and 180° respectively), where the following boundary 


conditions are satisfied. 
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At each plane of symmetry: 


ft- 0 

d((> 

ft-o 

(3.3) 

f =° 

dtp 

w = 0 

de 

d<p 


ae, 

-r— = o 

d(p 



At the cone surface 0 = d b : 


U = U C 

P = Pe 

(3.4) 

w — w e 

a P ^p e 
d<f> d<p 


o 

II 

> 

o 

II 

*1* 


At the main shock the conditions are given in terms of the free stream 

conditions and the 

ae s 

shock slope, ^ , 

<J<P 

by the Rankine-Hugoniot relations [8]. 



3.2. Viscous Flow 
3.2.1. Equations of Motion 

The viscous flow region, prior to separation from the body, is assumed to extend over 
a small distance normal to the cone body. The equations of motion for viscous, laminar, 
compressible, and nonisentropic flow about an inclined, axisymmetric, three dimensional 
body are nondimensionalized. The independent and dependent variables are written in terms 
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of the orthogonal coordinate system shown in (Fig. la.b). The Prandtl boundary layer 
assumptions are then applied to give the resulting equations of motion [18]. The Prandtl 


number (Pr) is assumed constant as are the specific heats c v and c p . 


Conservation of mass: 


a 

dx 





= 0 


(3.5a) 


Conservation of x momentum: 


3u du pw 3u 

pu S +pv 37 + "Ta^ 



dr 

dx 


dp c d du^ 

ax + dy ay J 


(3.5b) 


Conservation of <p momentum; 


dw dw pw dw 

pU 3r +pV 37 + T8» 


UW _ 1 dPc 

r dx r a^ 


3y 


' a w A 

p a y 


(3.5c) 


Conservation of energy: 


aH aH pw aH 1 a 

pu— +pv— + r - p r dy 



+ 




r 


\ 

M 

Pr 

d_ 

a 

^u" + w'j 


dy 

ay 

2 J 

> 


(3.5d) 


The velocities are nondimensionalized in the same manner as that used in the inviscid equa- 


tions of motion. Additional nondimensional variables are defined in the following manner. 
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r P°° 

r = F p°° - 

L Poo a 


x 

L 


U 0 


X — . Uoo s 


y = 


H = 


c p T + Vi | u 2 + w 5 


»— V 

poca L 


(3.6) 


where H is the total (stagnation) enthalpy. 

Equations (3.5a-d) can be simplified further. These governing equations of motion are 
three dimensional and include variable density, an explicit dependence on the local radius r, 
and severe gradients normal to the wall. By various transformations these effects can be 
eliminated or minimized in the equations [18], [48], [49]. 

It is assumed that the variation of viscosity across the boundary layer can be 
represented by a linear relation. The viscosity p is replaced as follows. 


fj. oo M oo 1 oo 

where C is some constant value. A full explanation of the equation is given by Chapman 
and Rubesin [50]. 

To remove the explicit dependence of the density we employ the Howarth transforma- 


tion. 
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x, = x. 



( 3 . 8 ) 


To reduce the effect of body curvature (r), the Mangier transformation is used. 


x 2 = Jrdx,, 
o 


y: = r yi 


( 3 . 9 ) 


This transformation also relates axially symmetric flow to plane flow. 


The transformed equations, at this point, resemble the equations of motion for 


incompressible flow in three dimensional Cartesian coordinates. 


We will further exploit the character of the equations in their simplified form. The 
outer edge of the boundary layer grows in a parabolic manner in the x direction ( y: = x 2 J )• 
A similarity variable can therefore be constructed. 


n = 


U e 




[ 2^00 ) 


Vi 


Xi 


( 3 . 10 ) 


The Blasius-type of transformation reduces the number of independent variables from three 
to two. 


The equations of motion are simplified further still. The dependent variables u, w, H, 
and v are replaced by u 2 , w 2 , s, and v 2 , respectively. They are defined as follows. 


u 

u-> = — 

' u e 


( 3 . 11 ) 
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w 

Wi = 

■ u e 


s = 1 - 


Vi = 


2H(y- 1) 
(7+ 1) 


R 

2 


r 3Re r poo 1 

1/2 ( 

k 2p e sin0 b j 

\ 

w 2 

3p< 


pv + u 


Ut + 


p e sin0 b d<p 



U^r 

where Re r = 

(U oo 


The results of our various transformations and substitutions are equations (3.12a-d). 


Conservation of mass: 


I 3wt 3vt 

— - + = 1 5 W 2 ~ 1.5U 2 

sin0 b dip arj 


(3.12a) 


Conservation of x momentum: 


W-> 3ui 3 ui 3“u 2 t 

"+v->^r r = l2 7-I1W2U2 + W5 


sin6> b 3 <p *3 r) *3r] 2 


(3.12b) 


Conservation of <f> momentum: 


a 2 

~ w 2 2 

-^+V,^-^ = l 11 (l-S) + l«U2 + l 8 W2 + l2-— T-'lW2- w 2 u 2 

sin3 b 3 <p " 3rj 3n" 


(3.12c) 


Conservation of energy: 
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w 2 8s 9s . 3 : s 

— + v 2 — *3 -> 

sin# b o<t> dr) dr)' 


U 9 : 




9 2 


, t( u 2 > + yfrW 

9tt - 9n 


(3.1 2d) 


The simplification of the equations of motion is now complete. The four partial 
differential equations describing three-dimensional, viscous, compressible boundary layer, 
with explicit dependence on curvature, has been converted into a system of equations 
resembling the boundary-layer equations for two-dimensional incompressible flow in Carte- 
sian coordinates. 


The coefficients l b ... , 1 12 are functions of the freestream conditions and the circum- 
ferential position <t>. The coefficients are defined below. 


1 _ «[e 

u e sin0 b 9 <p u e 


(3.13) 


!-> = 


1.5C 


(y-i) 

(y+D 


+ 


(y+ 1)M; 


Vi 




(y-D 
(y+ 1) 


r \ 



Pr 


V ) 


h 


l 


5 


l 3p e ^ 

2p e sin9 b 9 0 


We 

U e 


\ 


J 


(y — 1) 3p e 
yp e sin0 b 9 <p 



21 


1 foe | 2 W c 1 

2yp e sin0 b d0 u e . 


(y— 1 ) foe 
2yp e sin0 b dip 


(y+1) dfoe 
2yu*p e sin0 b fo 


(y- 1) ^ 2 Pe 
2yp e sin0 b dip 2 


( y + 1)>6 

li i - 

2 (y- l)u‘ 

I 1 ^ W e 

12 u e sin0 b dip 

Using the results from the previous chapter on the method of weighted residuals we 
will transform the equations 3.12a-d from their differential form to an integral form. The 
Crocco transformation (3.14) will be applied to the integral equations to change the 
independent variables from 0 and q to 0 and U 2 - 



o 


Since the variable v 2 has no immediate physical significance, we seek weighted combina- 
tions of the equations to eliminate the explicit dependence on the variable. 


We assume that the general weighting function, f(u 2 ), vanishes at the outer edge of 
boundary layer ( q = oo or u 2 = 1 ). By multiplying (3.12a) by f(u 2 ) and adding to 

x(3.12b) and integrating the transformed equation from zero to one, with respect to 
du 2 



~>2 


U 2 , we obtain (3.15a). 



w 

0 


Wt 


( I 


■j. Wl J- 

du 2 = sin0 b | I5 j f — -dui - 1.5 J f du 2 
o T o T 

!■ df 


, u 2 


df 3 t 
du-> du-> 


dui - 1 1 j 


df w 2 

U-> 

duj " x 


du 2 + 1 


du! r 


dui 


\ 


) 


df(u 2 ) 

The integral of the sum of w 2 x(3.12a), w 2 x — : - x(3.12b), and f(u 2 )x(3.12c) 

(3.15b). 


d !■ w 2 

— [ f — du-> = sind b 
d * o T 


l \ w 2 

1 7 J f — -du 2 - 2.5 j fu 2 — du 2 

r\ T n " 


U" 


1 

f 

df w 2 . 

duo + 

— du 

j 

du 2 t 

0 


i,j ^ 

3t 

W-> — — du-> 

'1 du: 

■ du-> 


(1-s) 


\ , 3r ^ w 2 


j. 3 ‘w 2 
+ U f fr— — du 2 

n 3"Ui 


The integral of the sum of sx(3.12a), sx 


df(u->) 

— x(3.12b), and f(u 2 )x(3.l2d) 

du 2 


(3.15c). 


_ 3 _ 

dip 




(3.15a) 


is equal to 


(3.15b) 


is equal to 


(3.15c) 
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df 

du 2 


dr 

3u 2 


du 2 


3w-> 


du 


-(ru->)du 2 + 1 4 f f— — ( tw 2 “t ~ )du 2 

* (7U ) ”H ^ 


du-> 


\ 


j 


It is assumed that there is no surface injection. Therefore v, and hence v 2 , is zero at the 
wall. Since v 2 and f(u 2 ) are zero at the wall and the outer boundary layer edge, respectively, 
v 2 no longer appears explicitly in the equations. 


Recall from chapter 2 that if linear interpolation is to be used in the finite element 
method, the derivatives within the integral must be less than second order. The application 
of integration by parts to (3.15a-c), and dropping the numbered subscripts, results in 
(3.16a-c); the Dorodnitsyn [29] boundary layer formulation for the inclined cone. 



(3.16a) 



(3.16b) 
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\ 


°J 


A f f^du = sine b f 1 5 f f^du - 1.5 } fu-du 
B<P 0 T 0 T o T 


!• df B s J , r ai sw r ui 
— (l-> + h) I — t — du — li I u — du + ■" 

12 • J 0 du du 1J 0 du r q du t 


df sw . r df sw 2 


du 


[ df urdu-l 4 [ TW lr du + l2 i dT Ai7^ ST ^ du 
. du 4 J du du J 0 du du 


- * 3 f * 


ds 

du 


(3.16c) 


The four partial differential equations describing three-dimensional, viscous, compres- 
sible boundary layer, with explicit dependence on curvature, (3.5a-d) have been converted 
into a system of three integral equations. The Dorodnitsyn formulation of the equations of 
motion offers some significant advantages over even the simplifications in equations 
(3.12a-d). 

To start, the weak form of the equations are being solved. This should decrease the 
errors in calculation since they are being averaged out across the boundary layer. By using 
u as an independent variable across the boundary layer, the infinite domain of q has been 
replaced by a finite domain. The high gradient of the variables in y have been decreased 
significantly, which results in high resolution of the dependent variables near the wall. This 
is of particular importantance in turbulent flow. The use of u as an independent variable 
also allows a uniform grid to be employed that automatically follows the growth of the 
boundary layer. The variable v 2 does not appear explicitly in the equations and can be 
recovered later. The shear stress, r, along the cone generator, is solved for directly and 
should be particularly accurate. 
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3.2.2. Initial Conditions 

In order to solve the integral equation, initial values of the dependent variables must 
be known. These values can be determined at the windward ray of symmetry. 

At <p = 0 : 


w = 0 




(3.17) 


Again neglecting the numbered subscripts, equations (3.16a) and (3.16c) reduce to the fol- 
lowing algebraic relations. 


f r ! d U = sin8 b ( - 1 .5 ) nfidu + ], } -f 

J T 90 [ l * l du au J 


r df dt . ) 


(3.18a) 


f f 1 |^du = sin0 b - 1.5 j fu^du - (1 2 + 1,) ] 
i r d<t> ( } 0 r J du du 


df 3s 


du 


(3.18c) 


- U f —urdu + 1*> f — ^-^-(sr)du - I3 

H J du “ l du du 


df 3 


du 3u 


fr 


3s 

3u 


Since the cross flow velocity w is zero, (3.16b) must be differentiated with respect to <p. 


This reduces to (3.18b) at the windward ray. 


1 

J f 


1 

T 



sin0 b 

2 


r 1 dw \ df 3 , 3w , 

- 2 - 5 l uf 7^ du+l! l^sr ( ^ t)du 


du 3u dtp 


(3.18b) 
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df d . 9w 

^ t sr ( '^ )du+ 




du - 19 



+ '4 f T d “- 


1 




3.2.3. Boundary Conditions 

For the equations of motion (0^0) the following Dirichlet conditions hold. 
At the cone wall: 


w = 0 
u = 0 



At the outer edge of the boundary layer: 

w e 

w= — 

Uc 

u= 1 
s = 0 

At <p =0 the same Dirichlet type boundary conditions of (3.16a-c) hold, with the following 
additions. 

At the cone wall: 

dtp 

At the B.L. edge: 

3w _ 1 9w c 
dtp u e dtp 
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4. APPLICATION OF FINITE ELEMENTS 


4.1. Equations of Motion 

As was explained in chapter 2, to determine the nodal values in the finite element 
method, the general weighting function, f(u), must be replaced by a set of linearly indepen- 
dent functions f k (u) ( for k= 1 N ). 

For this problem the set f k (u) is written as: 

f k (u) = (l -u) 2 N k (u), ( 4 -0 

where N k (u) is the linear interpolation function, at a particular node k. The term (1 -u) is 
introduced to satisfy the requirement that f k (u) equal zero at the outer edge of the boundary 
layer. The squaring of (1 -u) prevents singularities from appearing in the integrals when the 
dependent variables are given their finite element interpretation. 

The trial solutions for the dependent variables w, r , and s are shown. The group finite 
element formulation, described in reference [51], is used in making the trial solutions of the 
various combinations of dependent variables. 


w 

r 


(1-u) 


I N,blj 


1 

r 


(l-u) 


I NM 


sw 

r 


(l-u) 


X NjMj 


w = uX Njb4j 

j 


(4.2) 
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r = (l -u)£ Njb5j 

J 

s = X N j b6 j 


T 


1 

(1 - U) 


I N J b7 J 



X 


(1 - U) 


I N , b8 , 


s^ 1_ 

x ~ ( 1 -U) * 


I N n 


wr = u( 1 - u)£ NjblOj 


SW* U“ 


T ( 1 -U) 


In, mi, 


st = (1 -u)£ Njbl2j 


The simultaneous imposition of the particular analytic variation of u within each element on 
the variables prevents the exact interrelationships from being satisfied except at the nodes or 
in the limit The use of u outside of the summation sign is to insure the correct 

behavior of the dependent variables without invalidating the interrelationships at the nodes 
on the boundaries. 

Substituting the representations of (4.2) into equations (3.16a-c) gives 


Scii, 


dbl, 

d(f> 


= sin0 b 


r i 5 £Cl kj bl J -1.5Xcl kj b7 J + 

j j 


(4.3a) 
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b X C2kjb5j - 1 1 ^C3 kj bl j + ]T | C3 kj b2 J 
j J i j 

for k = 1 , ... , N 


YC4 kj = sin0 b fl 7 IC4 kj b2 J -2.5XC4 kj bl J + 

j d * j j 

\ 

l 2 X C5 kj b 1 Oj - 1 , X C7 kj b2, + X C7 kj b8j + 
j j j 

1 , , XC8 kj b7j - I, , IC8 kj b9j + l 8 XC4 kj b7 J - 


21i ^ C10A k b5 k _,b4 k _, + ClOB k b5 k _,b4 k + C10C k b5 k b4 k _, + 
( C10D k + C10E k )b5 k b4 k + C10F k b5 k b4 k+1 + 


C10G k b5 k+1 b4 k + C10H k b5 k+1 b4 k+1 j j 
for k = 2 N 


, ... , -I 2 b5,b4, for k = 1 


(4.3b) 


Y C 1 k j— ^ = sin^b ( lsX cl kj b 3 j“ 1 - 5 S C 1 ki b 9 J _ 1 | X C 3 kj b3 J + 

i "♦ i J i 

£C3 kJ bl 1 , - 1 4 ICV5, + 1, J C2 k ,bl2, - 


(1 2 + 1 3 ) ^CllA k b5 k _,b6 k _, + CllB k b5 k _,b6 k + CllC k b5 k b6 k _,+ 
( Cl 1 D k + Cl lE k )b5 k b6 k + Cl !F k b5 k b6 k+) + Cl !G k b5 k+1 b6 k + 


(4.3c) 


Cl lH k b5 k+ |b6 k+ | j j -l 4 (ci2A k bl0 k _ 1 b4 k . 1 +C12B k bl0 k _,b4 k + C12C k bl0 k b4 k _ l + 
( C12D k + C12E k )blO k b4 k + C12F k blO k b4 k+) + 020,^10^^ + 


C12H k blO k+ |b4 k+) j j 
for k = 2 N 


( 




, - 1 _ 3 b5 1 


(b6 2 -b6|) 

Au 


for k = 1 , 
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where Cl kj through C12H k are found in appendix A. With the use of finite elements, the 
integral equations of (3.16a-c) have been transformed into a system of ordinary differential 
equations. The solution of blj , b2 J( and b3j, can be solved by a variety of numerical 
methods. 


4.2 Initial Conditions 

The application of the finite element method to the integral equations governing the 
initial conditions (3.18a-c) is done in an identical manner. The trial solutions for the depen- 
dent variables are similar to those in equation (4.2). 


1 dw 
t dtp 


(1 -u) 


I N jQI, 


3w 


(1-u) 


X n jQ 2 j 


A dw 
r dtp 


(1-u) 


X N /»i 


dw 

d(p 


X N p4, 


r=(l-u>^N j Q5j 

J 


s = X n jQ 6 j 

J 


T 


1 

d-u) 


x w 


'3w'’ 

dtp 


(1-u) 


X w 


(4.4) 
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A _ 1 

T (1-U) 


I N,Q9, 


t|^ = u(l- U )£N J Q10 J 
OP i 


'dw' 2 


dip 


= Y nqh. 

(1-U)f J J 


ST = (1 -u)£ NjQ12j 


The Galerkin finite element formulation of equations (3.18a-c) is as follows: 


- £cl kj Qlj + sine b 


for k = 1 N 


-1.5XCl kj Q7 j + l 2 XC2 kj Q5 j 


= 0 


-XC4 kj Q2 r 


sin0 h 


-2.5X C4 k J Q 1 j + l 2lC5 kj Q10 J - 


l,IC7 kj b2 J+ XC7 kj b8 J + l^CSkjQ?, 

J J J 

-l 9 XC8 kj Q9j + 1,0204^ 

J J 

-21 2 [ciOA k Q5 k _iQ4 k _, +C10B k Q5 k _,Q4 k + C10C k Q5 k Q4 k _, + 
( C10D k + C10E k )Q5 k Q4 k + C10F k Q5 k Q4 k+ , + 

C10G k Q5 k+1 Q4 k + C10H k Q5 k+1 Q4 k+1 j ) =0 
for k = 2 N 


(4.5a) 


(4.5b) 


(, . -1 2 Q5iQ 4,) =0 for k = 1 
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-Xcl kj Q3j + sine b 


' - 1.5 ICl^j - l 4 IC9 kj Q5 j + l 2 IC2 kj Q12 j - 
j i i 


(1 2 + 1 3 ) J Cl 1 A k Q5 k _,Q6 k _ , + Cl lB k Q5 k _iQ6 k + Cl lC k Q5 k Q6 k _, + 
( Cl 1 D k + C 1 1 E k )Q5 k Q6 k + C 1 1 F k Q5 k Q6 k+ , + C 1 lG k Q5 k+ ,Q6 k + 


Cl lH k Q5 k+ ]Q6 k+ i j j =0 
for k = 2, ... , N 


(4.5c) 


,-l 3 Q5 


(Q6,-Q6,) 

Au 


= 0 for k = 1 , 


where Cl kj through Cl lH k are the same used in the equations of motion (4.3a-c). 

Notice that (4.5a-c) are essentially nonlinear algebraic equations in Qlj, Q2 j( and Q3j. 
The equations can be solved by an assortment of iterative methods. 


4.3 Boundary Conditions 

At the wall (u = 0): 


w = 0 


s s w 

For <p = 0 we have in addition: 



Which give the following nodal values: 


b3| = s w bl 1 
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Q3 ! = s w Qi i 

Notice that since t gives a Neumann and not a Dirichlet type of boundary condition, 
the value of t at the surface is not required by either equations (4.3a-c) or (4.5a-c). 

At the outer edge of the boundary layer (u = 1). 
w e 

w = — 

Ue 


r = 0 
s = 0 

For cp = 0, we have in addition: 

3w _ We 1 3w e 

d<t> u c ' u e d<p 

which give the following nodal values: 
w c 

b2 N = bl n 


b3 N = 0 


1 3w, 


Q3 n = 0 

Again, please note that the value of r at the boundary layer edge need not have been 
specified. However, if (1 -u) where not used in the representation of r, the nodal value of 
b5fl would have approached zero as the number of nodes increased ( N— »°° ). This would 
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mean that b2 N and bl N would have had to satisfy two conflicting relationships at u = 1 ( 


b2 N w e 

= — and 

bl N - u e 


b2 N 

bl N 2 


-»0 ). 



5. SOLUTION OF EQUATIONS OF MOTION 


In this section the steps taken to solve the governing equations numerically, are 
shown. 

5.1 Initial Conditions 

The set of ordinary differential equations of (4.3a-c) can be solved by a marching 
algorithm. However, any marching scheme considered will require the values of the nodes ( 
bl r b2j, b3j ) at some initial circumferential position, £, to start. The first step in solving for 
these nodal values, is to solve the equations of motion at <p =0. The initial condition equa- 
tions (4.5a-c), as mentioned earlier, are a system of 3N-3 nonlinear algebraic equations that 
must be solved iteratively to give the initial nodal values. The nonlinear equations are 
solved numerically by a quasi-Newton iterative scheme. This scheme requires a reasonable 
guess of the 3N-3 initial values to begin the process. This becomes rather difficult if an 
appreciable number of elements are to be considered. 

An efficient method for the solution of the initial conditions is presented. Using a 
similarity transformation, the partial differential equations (3.12a-d) are reduced to a system 
of ordinary differential equations. 


^ =-^(Tf+ 1 >2g) — 


dn 

d 3 g i 

dn 


h 2 


dn' 


l 9 Sy 


It 

-A + 


3 , , . d-g df dg 

— f + lpg) , - ~T~ j ”M2 

2 dn 2 dT l d n 

1 9 ( 1 s w ) 


fitl 

1 

ll ° I 

df 


l 12 sin0 b 



l 12 sin0 b li 2 sin0 b 


(B.4a) 


(B.4b) 


d 2 A Pr 


dn' 


s w l 2 


3 . . , dA U d 2 f 


dn 2 dT i 


B.4c) 


where 
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df 




iisinO b 


dg 

dn 


s = s w ( 1 - A) 


v = -(yf+l | 2 g ) 


(B.2) 


The boundary conditions are: 
At r\ =0: 


jdf 

dn 


= 0 , 



A = 0 


At rj = °°: 


f = 0 


g = ° 


dn 

A = 1 


^=1 

dn 


The details of the transformation, which is similar to that used by Moore [18], [48] and 
Reshotko [49], are discussed in appendix B, 


Equations (B.4a-c) represent a relatively straight forward two point boundary value 


d i d 2 

problem. The solution requires that the correct values for -, — — ^ 

dn" dn 


, and be found at 
dn 


the wall which will cause — , and A to approach the value of 1 as n approaches 

dn dn 


infinity. The numerical scheme used in integration was a fourth order Runge-Kutta method 
with a step size of 0.01. Integration up to an n of 6, was found to be sufficient. 
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3w 

The dependent variables t, - — , and s are found from (B.2) at the values of u required 

op 

by (4.5a-c). The nodal values of Qlj ,Q2j, and Q3j are then determined from the finite ele- 
ment representations of the three variables (4.4). At the boundaries however, the nodes are 
determined as follows: 

At u = 0 ( q =0 ): 


Ql,= 


d 2 e 

l| 2 sin0 b — 
dr? - 


d-f 

{^ 2 ) 


Q2,= 


d'g 

l 12 sin0 b -— 

drr 

w 


Q3| = 


d~g 

s w l, 2 stne b — *- 

djT 


'd 2 f^ 


dr? 2 J 


At u = i ( n = 00 ) ; 


Q1 n - -1: 


l (2 sin0 b 

(yf+ li2g) 


1 i 2 sin0 b 
Q^n = ~h 7 

(~f+hzS) 
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Q3 n = 0 


These values are then used as the initial guesses for the quasi-Newton iteration. Con- 
vergence of the iteration method is quite rapid. As the number of nodes increases for the 
finite element formulation, the values at the nodes approach those values determined trom 
the similarity equations. For an appreciable number of nodes ( N>5 ) the iteration process 
can be skipped all together and the values determined from the similarity equations alone. 

The starting values for the equations of motion are determined as follows. 

At some small angle e measured from the windward ray: 

nW 

w(e) = e-^-(0) 

r(e) = t(0) 
s(e) = s(0) 

therefore 


blj = eQlj (51) 

b2j = (e) 2 Q2j 
b3, = eQ3j 

The value of e to be used, must be small enough so that the approximations made for 
t and s are valid. However e should be large enough that the values of b2j do not approach 
the round off error of the machine. The value found to work well in this study was an 


£ Of 1°. 
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5.2 Equations of Motion 

The marching algorithm to be used in the solution follows the theta method, outlined 
in references [45] and [46], and discussed in chapter 2. The ordinary differential equations 
of (4.3a-c) are replaced by an implicit finite difference approximation in <f>. 


Equations (4.3a-c) are rewritten as: 


fci kj Abi; +1 =A0 


/3Dl k +l +(1 -/3)Dl k J 


(5.2a) 


X, C4 kj Ab2j n+ 1 = A<p 


/3D2 k +1 + (1 - 



(5.2b) 


fci kj Ab3j n+l = A^ 

J 


/3D3 k +1 + (1 -/3)D3 k 


(5.2c) 


for k = 1 N where 


Ablj n+I = bl" +l - bl" 

Ab2j n+1 = b2j n+l -b2j n 
Ab3j n+I = b3" +1 - b3" 

Dl k , D2 k , and D3 k are the right hand sides of equation (4.3a), (4.3b) and (4.3c), respec- 
tively. The superscript n denotes a particular circumferential position. Here, the parameter 

P (0^/3 <1) is introduced to control the degree of implicitness. The quanities Dl k +I , D2 k +I , 

and D3 k +I are linearized by expansion about the known position <p n [52]. 
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Dl k +1 = Dl k (0 n+1 ) + 
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D3 k +I = D3 k (0 n+I ) + 
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The linearization is substituted into (5.2a-c) and results in the following system. 
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(5.4b) 


(5.4c) 


DDl k , DD2 k , and DD3 k are Dl k , D2 k , and D3 k respectively using the modified coefficients 

lXj, ... , lX|2- 


For example 
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lx, = /3lT +1 + (1-0)1" 
dDl k ) " f9Dl k )" 

The arrays — , -r are Jacobians. The integrals Cl kj C12H k are 

obi, ob2, 

V i) V J 

functions only of u and therefore need only be evaluated once. The integrands reduce to 
polynomials of various orders and can be solved to very high accuracy by Gaussian quadra- 
ture [44], The highest order polynomial encountered is six and therefore quadrature of an 
order of no higher than three need be used to solve all of the integrals. 


The equations (5.4a-c) can be lumped into a single matrix equation as shown. 

3N 

£LHS im Ab m n+, = RHS, 1=1 3N (5.5) 

m=l 

where 

Ab m n+l = Abl j n+, ,Ab2 j n+, ,Ab3 j n+1 j=l, ... ,N 

b m = bl Jt b2j,b3j 

The array LHS lm is a sparse 3Nx3N matrix and RHS, is a 3N vector containing the vectors 
A</>DDl k , A<f>DD2 k , and Atf>DD3 k respectively. 


The boundary conditions are as follows: 

Ab 2N n+1 =l, n+l Ab N n+l + (l, n+l -l,j bN n 
Ab 2 N+1 n+l = s w Ab N n+1 


(5.6) 
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By incorporating the boundary conditions into equation (5.5) through the use of Lagrange 
multipliers [27], the length of the vector RHS, is increased by three. The array LHS )m is 
now a (3N+3)x(3N+3) matrix. 

Equation (5.5) is solved at each step using a general LU decomposition algorithm. 
Rather than iterate at a particular <p position, the method varies the step size A <p to achieve 
the desired accuracy. Before computing the new solution bj n+l = Abj n+I + bj" the step-size is 

Ab, n+I 

determined in the following manner. As long as the maximum value of the ratio, — - — , 

bj" 

exceeds a preset constant the step-size is halved, if this step-size is greater than the preset 
minimum. If the maximum value of the ratio is less than one tenth of the preset constant, 
the step-size is increased by 50%, if this step-size is less than the preset maximum. Since 
no iteration is required to calculate the solution at each <p, and since only moderate effort is 
required for the solution of the implicit difference equation, considerable computational 
efficiency is achieved. 
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6. RESULTS 


6.1 Convergence Properties 

The convergence rales for finite element approximations of w, t, and s are plotted in 
figures 4, 5, and 6. The ordinate is the measure of the average discrete L2 norm of the rela- 
tive error. The error used is calculated in the following manner. 
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error(s) = 
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U= 2 

v S J" J 

l 



(N-2) 


(6.1c) 


Where w jex , r jex , and s jex are the exact values at the respective node. 

In order to eliminate errors due to the discretization of 0, in the spatial convergence 
results, the convergence study was taken at the windward ray. Recall from chapter 5 and 
appendix B, that the windward equations can be represented by the three ordinary 
differential equations 


d*f 
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A solution to these equations, for a = 10°, were obtained using a fourth order Runge- 
Kutta scheme with An =0.01. This solution is used as the exact values in the convergence 
results. 

The theoretical global convergence rate for linear elements applied to a linear problem 
is first order. It was expected that since the cone problem involved highly nonlinear equa- 
tions and the solution involved the simultaneous prescription of w, r, and s at each node, 
that the convergence rate would be less than first order. This is seen for the variables s and 
r. However, the variable w achieves almost second order convergence. The errors of w, r, 
and s for crude grids are quite good. For example with only two elements s has an error of 
approximately 15%, w has 11% and x has only 5%. For four elements the errors for s, w, 
and r are 1 1%, 2%, and 3% respectively and for eight elements s.w, and t are 6%, 1%, and 
2% respectively. Therefore, even though the theoretically expected convergence rate is not 
achieved for s and r, this is offset by the high accuracy on course grids. For the variable w 
the benefits of both high convergence and accuracy on course grids, are enjoyed. 


6.2 Numerical Results 

The results presented are for a sharp cone with a half angle of 15° at angles of attack 
of 5° and 10°. The fluid properties used are for air. The flow parameter are given below. 
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0 b = 15° (6-2) 
Re x = 5xl0 5 
M« = 7 
Pr = 1 


C = 0.922 


Po c = 0.282 

U oe = 2.333 



Only two elements ( Au = 0.5 ) were used in the calculations. Through numerical experi- 
mentation the value of 10~ 2 was used to ensure stability during <p marching. The step size ( 
A </> ) varied from 10~ 5 to 10 _: \ This program was run in double precision and on a VAX 
11/780 at the NASA Ames Research Center. 

Figure 7 shows the distribution of the surface coefficient of friction ( C fw ) with the 
circumferential position <p for a = 5° and 10° respectively. The coefficient is defined below. 


Qw — 


_ 3w 

Kp ooUo 


u 0 


6p«, 

R e xPc 


a w 

1 T a^T 




(6.3) 


The influence of incidence on C fw can be seen by comparing the curves of figure 7. For 
both angles of attack the C fw varies smoothly from 0=0° reaching a maximum value at 
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approximately <p = 100° and falling off again toward zero as the leeward ray is approached. 
The values of Cf w for a = 10° are greater than that for a = 5° at all circumferential posi- 
tions. The failure of C fw to equal zero at <f> = 180° for both angles of attack is probably 
caused by the buildup of interpolation errors, from the crude grid used (two elements), dur- 
ing marching. 

Figure 8 shows the variation of the surface coefficient of friction in the radial direction 
( C fu ) with <p. The coefficient Cf u is defined as follows. 


Q u - 


r)u 

h 


'/2 Poo U c 



(6.4) 


As can be seen in the figure the maximum value of C fu occurs at the windward ray ( </> = 0° 
) and falls smoothly to a nonzero value as <f> approaches the leeward ray. The distribution of 
C fu is dependent on how the velocity u e and the thickness of the boundary layer varies cir- 
cumferentially. For the test cases u e increases with </> which would increase Cf u . However, 
the increase in boundary layer thickness, which decreases Cf u , overcomes the influence of 
u e . The angle a =10° is seen to give greater values of C fu than a = 5° at <f> =0°. This 
difference in values of C fu decreases as <p advances. The decrease in difference seems to be 
a result of the insensitivity of the external conditions to incidence except in the windward 
region. 


The relative heat transfer at the wall 


Qw 


) is plotted, versus <f>, in figure 9. The 


heat transfer is calculated in the following way. 
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The relative heal transfer seems to follow the same trends as the coefficient C fu . As was 


seen for C, u the maximum value of 



occurs at the windward ray and decreases 


smoothly as <p increases. Also, as was seen for C fu , the effect of incidence is evident mostly 
in the windward region. These similar trends between heat transfer and C fu were expected 
since the Prandtl number (Pr) was set equal to 1. The flow variable Pr = 1 makes the 
momentum and thermal boundary layers identical. 


Figures 10 and 11 show the cross flow velocity profiles within the boundary layer at 
0 =45°, 90°, and 135°, for a = 5° and 10° respectively. The growth of the boundary layer 
with increasing 0 is apparent. The profile curves are typical of three dimensional boundary 
layer flows, with a cross flow velocity maximum occuring within the layer [21]. Also, the 
influence of incidence on the profiles is noticable. As a increases from 5° to 10° the cross 
flow velocity increases, the velocity maximum shown at 0 = 135° becomes more exag- 
gerated, and the boundary layer increases in thickness. 

The effects of incidence given are similar to the effects shown in references 


[l],[20],[21j,[53]. 
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7. CONCLUSIONS 


7.1 Summary 

A numerical method has been proposed to model the three-dimensional boundary 
layer about a yawed circular cone in a supersonic stream. It can, in principal, be extended 
to apply to the boundary layer of pointed bodies of arbitrary cross section. 

The method reduces the three-dimensional boundary layer equations to two- 
dimensional form, through a coordinate transformation, and eliminates the explicit appear- 
ance of the normal velocity. The coordinate transformation also allows high resolution in 
the viscous region near the cone surface. The reduced boundary layer formulation is then 
presented in a Petrov-Galerkin finite element form and discretized across the layer using 
linear interpolation functions. Linear interpolation is the most computationally efficient to 
use from the family of piecewise continuous polynomials. The finite elements yield a sys- 
tem of ordinary differential equations in the circumferential coordinate. The system of 
differential equations, when put in a matrix form, yield a sparse mass matrix. The circum- 
ferential derivatives are then solved by a noniterative implicit marching scheme that gives 
both speed and stability. The results shown are in keeping with those given in previous stu- 
dies. The proposed method gives acceptable accuracy with a very crude grid. 

7.2 Future Work 

Future applications of the finite element / finite difference method, by the author, will 
probably proceed as follows. 

(i) The semidiscrete Galerkin modelling of the boundary layer about the cone will be 
linked with an interacting numerical model of the inviscid conical flow [24] for high 
angle of attack studies. 

(ii) In a manner similar to that given by Holt [54]-[57], the boundary layer modelling will 
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also be carried on beyond flow separation into the reversed flow region and its reat- 
tachment. 

(iii) The complete method will be applied to various cross sectional geometries. 
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Fig. 1. Coordinate system and velocity components. 
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APPENDIX A 


A Definition of Cl kj Through C12H k 

The components Cl kj through C\2H± first mentioned in chapter 3, are defined below. 
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(1-u) -2N k 

du 


U k 


dN k 

du 


N k du 


dN k 

du 


N k du 


dN k+ , 

du 


N k du 


dN k 

N k+ |du 

du k 


dN k+l 

du 


N k+ |du 


= | u(l-u) 2 


u»-i 


dN k 

(1-u) — “2N k 
du 


dN k-i KI 

u— +N k _, 

du 


N k _,du 


= J u(l-u ) 2 


u t -i 


dN k 

(1-u)— -2N k 
du 



dN k 


u + N k 


du 


N k _,du 


U k 


= j u(l-u ) 2 


= J u(l-u ) 2 


u»-l 


dN k 

(1-u)— 2N k 

du 


dN k 

(1-u)— --2N k 


du 



dN k _, 


u — — — + N 

- 

du 

-| 

dN k 


u + N k 


du 


k-1 


N k du 


N k du 


U**i 


= J u(l-u ) 2 


dN k 

(1-u)— -2N k 
du 



dN k 


u — + N k 


du 


N k du 


Ukfj 


= J u(l-u ) 2 


u» 


dN k 

(1-u)— -2N k 
du 


dN k+ , 

u— +N k+ , 

du 


N k du 
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V , dN k 

C12G k = J u(l-u)‘ ( 1 -u)— -2N k 

U4 L 

Uk + 1 dN k II dN k+) _ I , 

C12H k = u(l-u)- (1-u)- 2N k u— — + N k+ | N k+) du 

J du du 

Ut L J L J 


dN k 

u— — + N k N k+ |du 
du 


II 


7f 

If 




